#include "jrm.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#define IS_INITIALIZED  (MA != NULL)
#define WAS_EXECUTED    (VX != NULL)

#define ABS(X)          ((X) > 0 ? (X) : -(X))

/* Matrices and vectors order */
static int ORDER    = -1;
/* Row to be tested */
static int ROW_TEST = -1;
/* Iterations limit */
static int ITE_MAX  = -1;
/* Current iteration */
static int ITER     = -1;

/* Error precision */
static float ERROR  = 0;

/* Matrix A */
static float **MA = NULL;
/* Vector b */
static float  *VB = NULL;

/* Matrix A* */
static float **MA_S = NULL;
/* Vector b* */
static float  *VB_S = NULL;

/* Vector x */
static float *VX = NULL;

void jrm_initialize(char *input)
{
    int i, j;

    if (IS_INITIALIZED)
        jrm_delete();

    /* Parse input file */
    ORDER = (int) strtol(input, &input, 10);
    ROW_TEST = (int) strtol(input, &input, 10);
    ERROR = (float) strtod(input, &input);
    ITE_MAX = (int) strtol(input, &input, 10);

    assert(ORDER > 0);

    /* Allocate memory for A, A*, b and b* */
    MA = (float **) malloc(ORDER * sizeof(float *));
    MA_S = (float **) malloc(ORDER * sizeof(float *));

    VB = (float *) malloc(ORDER * sizeof(float));
    VB_S = (float *) malloc(ORDER * sizeof(float));

    /* Fill matrices and vectors */
    for (i = 0; i < ORDER; i++)
    {
        MA[i] = (float *) malloc(ORDER * sizeof(float));
        MA_S[i] = (float *) malloc(ORDER * sizeof(float));

        for (j = 0; j < ORDER; j++)
            MA[i][j] = (float) strtod(input, &input);
    }

    for (i = 0; i < ORDER; i++)
        VB[i] = (float) strtod(input, &input);
}

void jrm_delete()
{
    int i;

    if (!IS_INITIALIZED)
        return;

    /*
     * Free all memory and set original values in variables.
     */

    for (i = 0; i < ORDER; i++)
    {
        free(MA[i]);
        free(MA_S[i]);
    }

    free(MA);
    free(MA_S);
    free(VB);
    free(VB_S);

    if (WAS_EXECUTED)
        free(VX);

    ORDER    = -1;
    ROW_TEST = -1;
    ITE_MAX  = -1;
    ITER     = -1;

    ERROR  = 0;

    MA = NULL;
    VB = NULL;
    MA_S = NULL;
    VB_S = NULL;
    VX = NULL;
}

static void jrm_execute()
{
    int i, j;
    float x, error;

    /* Allocate memory for vector x */
    VX = (float *) malloc(ORDER * sizeof(float));

    /* Calculate matrix A* using matrix A */
    for (i = 0; i < ORDER; i++)
        for (j = 0; j < ORDER; j++)
            if (i != j)
                MA_S[i][j] = MA[i][j] / MA[i][i];
            else
                MA_S[i][j] = 0;

    /* Calculate vector b* using vector b and matrix A */
    for (i = 0; i < ORDER; i++)
    {
        VB_S[i] = VB[i] / MA[i][i];
        VX[i] = VB_S[i];
    }

    /* Iterate up to ITE_MAX */
    for (ITER = 0; ITER < ITE_MAX; ITER++)
    {
        error = 1e10;

        /* Update values of the vector x */
        for (i = 0; i < ORDER; i++)
        {
            x = VB_S[i];

            for (j = 0; j < ORDER; j++)
                x -= MA_S[i][j] * VX[j];
            x += MA_S[i][i] * x;

            /* Update error */
            if (ABS(x - VX[i]) / ABS(x) > error)
                error = ABS(x - VX[i]) / ABS(x);

            VX[i] = x;
        }

        /* Check stop condition */
        if (error < ERROR)
            break;
    }
}

void jrm_evaluate(int *iter, int *row_test, float *estimated, float *real)
{
    int i;

    if (!IS_INITIALIZED)
        return;

    if (!WAS_EXECUTED)
        jrm_execute();

    /* Calculate estimated value in ROW_TEST */
    *estimated = 0;
    for (i = 0; i < ORDER; i++)
        *estimated += MA[ROW_TEST][i] * VX[i];

    *real = VB[ROW_TEST];
    *row_test = ROW_TEST;
    *iter = ITER;
}
